#### setting environment ####
require(estimatr)
require(interflex)

#### data preparation ####
data <- read.csv("replication_data.csv")

## number of respondents
nrow(data)

## variables
# approval for Japan's participation in the TPNW
data$approve <- 6 - data$outcome
# dummy for support
data$support.bin <- (data$approve > 3) * 1
# dummy for opposition
data$oppose.bin <- (data$approve < 3) * 1
# security threats treatment
data$security <- (data$group == 2 | data$group == 4) * 1
# international pressure treatment
data$pressure <- (data$group == 3 | data$group == 4) * 1
# dummy for female respondents
data$female <- ifelse(data$gender == 3, NA, data$gender - 1)
# dummy for conservative respondents
data$conservative <- (data$self.placement > 6) * 1
# dummy for liberal respondents
data$liberal <- (data$self.placement < 6) * 1
# dummies for randomization blocks (used for interflex())
block.dummies <- diag(27)[data$block, ][, -1]
colnames(block.dummies) <- paste0("block.", 1:ncol(block.dummies))
data <- cbind(data, block.dummies)

#### census distribution (Online Appendix C) ####
census.age.gender <- read.csv("census_age_gender.csv")
census.prefecture.age <- read.csv("census_prefecture_age.csv")
census.age.education <- read.csv("census_age_education.csv", 
                                 na.strings = "-")

## Table A.1 (gender)
# survey
round(prop.table(table(data$gender)), 3)
# census
round(prop.table(colSums(census.age.gender[census.age.gender$age >= 18, 
                                           c("male", "female")])), 3)

## Table A.1 (age)
# survey
round(prop.table(c(sum(data$age < 30), 
                   sum(data$age >= 30 & data$age < 40), 
                   sum(data$age >= 40 & data$age < 50), 
                   sum(data$age >= 50 & data$age < 60), 
                   sum(data$age >= 60))), 3)
# census
round(prop.table(c(sum(census.age.gender[census.age.gender$age >= 18 & 
                                           census.age.gender$age < 30, -1]), 
                   sum(census.age.gender[census.age.gender$age >= 30 & 
                                           census.age.gender$age < 40, -1]), 
                   sum(census.age.gender[census.age.gender$age >= 40 & 
                                           census.age.gender$age < 50, -1]), 
                   sum(census.age.gender[census.age.gender$age >= 50 & 
                                           census.age.gender$age < 60, -1]), 
                   sum(census.age.gender[census.age.gender$age >= 60, -1]))), 3)

## Table A.1 (education)
# survey
round(prop.table(c(sum(data$education < 3), 
                   sum(data$education >= 3 & data$education < 6), 
                   sum(data$education >= 6))), 3)
# census
census.age.education[1, -1] <- census.age.education[1, -1] * 
  (sum(census.age.gender[census.age.gender$age >= 18 & 
                           census.age.gender$age < 20, -1]) / 
     sum(census.age.gender[census.age.gender$age < 20, -1]))
census.age.education[1, c("HS", "JC", "C", "GS")] <- 
  census.age.education[1, c("HS", "JC", "C", "GS")] + 
  census.age.education[1, "student"] * 
  (census.age.education[3, c("HS", "JC", "C", "GS")] / 
     sum(census.age.education[3, c("HS", "JC", "C", "GS")]))
census.age.education[2, c("JC", "C", "GS")] <- 
  census.age.education[2, c("JC", "C", "GS")] + 
  census.age.education[2, "student"] * 
  (census.age.education[3, c("JC", "C", "GS")] / 
     sum(census.age.education[3, c("JC", "C", "GS")]))
round(prop.table(c(sum(census.age.education[, c("ES", "JHS", "HS")]), 
                   sum(census.age.education[, c("JC")]), 
                   sum(census.age.education[, c("C", "GS")], na.rm = TRUE) + 
                     sum(census.age.education[3:11, c("student")]))), 3)

## Table A.1 (region)
# survey
round(prop.table(c(sum(data$prefecture %in% c(2, 3, 8, 12, 16, 24, 45)), 
                   sum(data$prefecture %in% c(4, 10, 14, 19, 35, 39, 41)), 
                   sum(data$prefecture %in% c(1, 6, 9, 15, 26, 29, 38, 43, 47)), 
                   sum(data$prefecture %in% c(13, 22, 23, 28, 33, 36, 44)), 
                   sum(data$prefecture %in% c(5, 11, 17, 20, 31, 37, 40, 42, 46)), 
                   sum(data$prefecture %in% c(7, 18, 21, 25, 27, 30, 32, 34)))), 3)
# census
round(prop.table(c(sum(census.prefecture.age[census.prefecture.age$prefecture %in% 
                                               c("Hokkaido", "Aomori", "Iwate", 
                                                 "Miyagi", "Akita", "Yamagata", 
                                                 "Fukushima"), -1]), 
                   sum(census.prefecture.age[census.prefecture.age$prefecture %in% 
                                               c("Ibaraki", "Tochigi", "Gunma", 
                                                 "Saitama", "Chiba", "Tokyo", 
                                                 "Kanagawa"), -1]), 
                   sum(census.prefecture.age[census.prefecture.age$prefecture %in% 
                                               c("Niigata", "Toyama", "Ishikawa", 
                                                 "Fukui", "Yamanashi", "Nagano", 
                                                 "Gifu", "Shizuoka", "Aichi"), -1]), 
                   sum(census.prefecture.age[census.prefecture.age$prefecture %in% 
                                               c("Mie", "Shiga", "Kyoto", 
                                                 "Osaka", "Hyogo", "Nara", 
                                                 "Wakayama"), -1]), 
                   sum(census.prefecture.age[census.prefecture.age$prefecture %in% 
                                               c("Tottori", "Shimane", "Okayama", 
                                                 "Hiroshima", "Yamaguchi", 
                                                 "Tokushima", "Kagawa", "Ehime", 
                                                 "Kochi"), -1]), 
                   sum(census.prefecture.age[census.prefecture.age$prefecture %in% 
                                               c("Fukuoka", "Saga", "Nagasaki", 
                                                 "Kumamoto", "Oita", "Miyazaki", 
                                                 "Kagoshima", "Okinawa"), -1]))), 3)

#### factual manipuration checks (FMCs) ####
## correct answer rate for the FMC for the security threats treatment
round(mean(data$MC.1 == 2, na.rm = TRUE), 3)

## correct answer rate for the FMC for the international pressure treatment
round(mean(data$MC.2 == 4, na.rm = TRUE), 3)

#### distribution of the outcome variable ####
## entire distribution
round(prop.table(table(data$approve)), 3)

## Figure 1(a)
png("Figure_1a.png", 
    width = 5, height = 3, units = "in", pointsize = 9, res = 1200)
par(mar = c(2, 3.5, 2, 0.5), lwd = 0.5)
plot(NULL, NULL, type = "n", bty = "n", xlim = c(0.5, 5.5), ylim = c(0, 0.5), 
     main = "Should Japan join the TPNW?", 
     xlab = "", ylab = "", xaxt = "n", yaxt = "n")
abline(h = seq(0, 0.4, 0.1), lty = 3, col = "gray80")
hist(data$approve, breaks = seq(0.5, 5.5, 1), freq = FALSE, 
     bty = "n", xlim = c(0.5, 5.5), ylim = c(0, 0.4), 
     main = "", xlab = "", ylab = "", xaxt = "n", yaxt = "n", add = TRUE)
segments(0.55, 0.44, 2.45, 0.44)
segments(0.55, 0.43, 0.55, 0.44)
segments(2.45, 0.43, 2.45, 0.44)
text(1.5, 0.44, 
     paste0("Con-TPNW\n", 
            sprintf("%4.1f", mean(data$approve < 3) * 100), "%"), 
     pos = 3)
segments(2.55, 0.44, 3.45, 0.44)
segments(2.55, 0.43, 2.55, 0.44)
segments(3.45, 0.43, 3.45, 0.44)
text(3, 0.44, 
     paste0("Undecided\n", 
            sprintf("%4.1f", mean(data$approve == 3) * 100), "%"), 
     pos = 3)
segments(3.55, 0.44, 5.45, 0.44)
segments(3.55, 0.43, 3.55, 0.44)
segments(5.45, 0.43, 5.45, 0.44)
text(4.5, 0.44, 
     paste0("Pro-TPNW\n", 
            sprintf("%4.1f", mean(data$approve > 3) * 100), "%"), 
     pos = 3)
axis(2, at = seq(0, 0.4, 0.1), labels = seq(0, 40, 10), lwd = 0.5)
mtext(c("No", "Probably not", "Can\U02BCt say", "Probably should", "Yes"), 
      1, 0, at = 1:5)
mtext("Percentage", 2, 2.5, at = 0.2)
dev.off()

## distribution in each group
# minimum pro-TPNW rate in experimental groups
round(min(tapply(data$approve > 3, data$group, mean)), 3)
# maximum con-TPNW rate in experimental groups
round(max(tapply(data$approve < 3, data$group, mean)), 3)

## Figure A1
# group means
group.mean <- tapply(data$approve, data$group, mean)

group.labels <- c("Group C (no information)", "Group T1 (security environment)", 
                  "Group T2 (international pressure)", "Group T3 (both information)")

cairo_pdf("Figure_A1.pdf", width = 6, height = 4, pointsize = 7)
layout(matrix(1:4, 2, 2, byrow = TRUE))
par(mar = c(3.5, 3.5, 3, 0.5), lwd = 0.5)
for (i in 1:4) {
  plot(NULL, NULL, type = "n", bty = "n", xlim = c(0.5, 5.5), ylim = c(0, 0.5), 
       main = group.labels[i], xlab = "", ylab = "", xaxt = "n", yaxt = "n")
  abline(h = seq(0, 0.4, 0.1), lty = 3, col = "gray80")
  hist(data$approve[data$group == i], breaks = seq(0.5, 5.5, 1), freq = FALSE, 
       bty = "n", xlim = c(0.5, 5.5), ylim = c(0, 0.4), 
       main = "", xlab = "", ylab = "", xaxt = "n", yaxt = "n", add = TRUE)
  segments(0.55, 0.44, 2.45, 0.44)
  segments(0.55, 0.43, 0.55, 0.44)
  segments(2.45, 0.43, 2.45, 0.44)
  text(1.5, 0.44, 
       paste0("Con-TPNW\n", 
              sprintf("%4.1f", mean(data$approve[data$group == i] < 3) * 100), "%"), 
       pos = 3)
  segments(2.55, 0.44, 3.45, 0.44)
  segments(2.55, 0.43, 2.55, 0.44)
  segments(3.45, 0.43, 3.45, 0.44)
  text(3, 0.44, 
       paste0("Undecided\n", 
              sprintf("%4.1f", mean(data$approve[data$group == i] == 3) * 100), "%"), 
       pos = 3)
  segments(3.55, 0.44, 5.45, 0.44)
  segments(3.55, 0.43, 3.55, 0.44)
  segments(5.45, 0.43, 5.45, 0.44)
  text(4.5, 0.44, 
       paste0("Pro-TPNW\n", 
              sprintf("%4.1f", mean(data$approve[data$group == i] > 3) * 100), "%"), 
       pos = 3)
  arrows(group.mean[i], 0.05, group.mean[i], 0.01, length = 0.03)
  text(group.mean[i], 0.04, 
       paste0("Avg\n", sprintf("%4.2f", group.mean[i])), pos = 3)
  axis(1, at = 1:5, FALSE, lwd = 0.5)
  axis(2, at = seq(0, 0.4, 0.1), labels = seq(0, 40, 10), lwd = 0.5)
  mtext(c("1\nShould not join\n", "2\nShould probably\nnot join", 
          "3\nI cannot say\neither way", "4\nShould probably\njoin", 
          "5\nShould join\n"), 1, 2.25, at = 1:5, cex = 0.65)
  mtext("Percentage", 2, 2.5, at = 0.2, cex = 0.8)
}
dev.off()

## relationship between the outcome and ideological self-placement (Figure 1(b))
mosaicplot.data <- table(data$self.placement, data$approve)
rownames(mosaicplot.data) <- 0:10
colnames(mosaicplot.data) <- rep("", 5)
cumsum(rev(prop.table(table(data$approve))))

png("Figure_1b.png",
    width = 6, height = 4, units = "in", pointsize = 10, res = 1200)
par(mar = c(0.5, 0, 0.5, 0), oma = c(0, 0, 1.5, 0), lwd = 0.5)
mosaicplot(mosaicplot.data, main = "", 
           color = c("gray50", "gray65", "gray80", "gray95", "white"), 
           border = "gray50")
text(rep(0.425, 5), c(0.12, 0.43, 0.725, 0.915, 0.97), 
     c("Should join", "Should probably join", "I cannot say either way", 
       "Should probably not join", "Should not join"), 
     cex = 0.9)
mtext("Ideological self-placement", side = 3, line = 0.5, 
      outer = TRUE)
mtext(c("Liberal (left)", "Conservative (right)"), side = 3, line = 0.25, 
      at = c(0, 1), adj = c(0, 1), cex = 0.8)
mtext("Approval of Japan\U02BCs accession to the TPNW", side = 2, line = -1, 
      outer = TRUE)
dev.off()

#### test of Hypotheses 1 and 2 ####
## estimation of Model (1) (Table 2, columns 1-3)
model.1 <- lm_robust(approve ~ security + pressure, 
                     data, fixed_effects = block)
round(summary(model.1)$coefficients[, c(1, 2, 4)], 3)

## Figure 2
group.CI <- tapply(data$approve, data$group, function(x) t.test(x)$conf.int)

png("Figure_2.png",
    width = 5, height = 3, units = "in", pointsize = 7, res = 1200)
layout(matrix(1:2, 1, 2), widths = c(6, 4))
par(mar = c(2, 3.5, 2, 2), lwd = 0.5)
plot(NULL, NULL, bty = "n", xlim = c(0.5, 4.5), ylim = c(3.5, 4.2), 
     main = "Mean approval of each group", 
     xlab = "", ylab = "", xaxt = "n", yaxt = "n")
abline(h = seq(3.5, 4.2, 0.1), lty = 3, col = "gray80")
for (i in 1:4) {
  segments(i, group.CI[[i]][1], i, group.CI[[i]][2])
  points(i, group.mean[i], pch = 19)
}
axis(2, lwd = 0.5)
mtext(c("Group C", "Group T1", "Group T2", "Group T3"), 1, 0, at = 1:4)
mtext("Approval of Japan's accession to the TPNW", 2, 2.5)
plot(NULL, NULL, bty = "n", xlim = c(0.5, 2.5), ylim = c(-0.4, 0.3), 
     main = "Average treatment effect", 
     xlab = "", ylab = "", xaxt = "n", yaxt = "n")
abline(h = 0, lty = 1, col = "gray80")
abline(h = seq(-0.4, 0.3, 0.1)[-5], lty = 3, col = "gray80")
segments(1, model.1$conf.low["security"], 
         1, model.1$conf.high["security"])
points(1, model.1$coefficients["security"], pch = 19)
segments(2, model.1$conf.low["pressure"], 
         2, model.1$conf.high["pressure"])
points(2, model.1$coefficients["pressure"], pch = 19)
axis(2, lwd = 0.5)
mtext(c("Security\nenvironment", "International\npressure"), 1, 1, at = 1:2)
mtext("Average treatment effect", 2, 2.5)
dev.off()

## substantive interpretation (Footnote 16)
# average outcome in the control group
round(mean(data$approve[data$group == 1]), 3)
# reduction rate in the extent of support
round(model.1$coefficients["security"] / 
        (mean(data$approve[data$group == 1]) - 1), 3)

#### test of Hypotheses 3 and 4 ####
## estimation of Model (2) with responents' gender (Table 2, columns 4-6)
model.2.gender <- lm_robust(approve ~ security + pressure + 
                              (security + pressure) : female, 
                            data, fixed_effects = block)
round(summary(model.2.gender)$coefficients[, c(1, 2, 4)], 3)
# The number of observations is less than in other analyses because 
# we excluded respondents who chose the “neither” option for their gender.
model.2.gender$nobs

## estimation of Model (2) with responents' ideological self-placement (Table 2, columns 7-9)
model.2.ideology <- lm_robust(approve ~ (security + pressure) * self.placement, 
                              data, fixed_effects = block)
round(summary(model.2.ideology)$coefficients[, c(1, 2, 4)], 3)

## treatment effects of the security threats treatment conditioned by ideological self-placement
# conditional average treatment effects (CATEs)
CATE.est.ideology <- model.2.ideology$coefficients["security"] + 
  c(1:11) * model.2.ideology$coefficients["security:self.placement"]
round(CATE.est.ideology, 3)

# reduction rate in the extent of support for each category of ideological self-placement
round(CATE.est.ideology / (mean(data$approve[data$group == 1]) - 1), 3)

## Figure 3
# standard errors of the CATEs
CATE.SE.ideology <- sqrt(model.2.ideology$vcov["security", "security"] + 
                           model.2.ideology$vcov["security:self.placement", "security:self.placement"] * c(1:11) ^ 2 + 
                           2 * model.2.ideology$vcov["security", "security:self.placement"] * c(1:11))
# distribution of responents' ideological self-placement
hist.self.placement <- prop.table(table(data$self.placement))

png("Figure_3.png",
    width = 5, height = 3, units = "in", pointsize = 7, res = 1200)
par(mar = c(3.5, 3.5, 0.5, 0.5), lwd = 0.5)
plot(NULL, NULL, bty = "n", xlim = c(0.5, 11.5), ylim = c(-0.6, 0.4), 
     xlab = "", ylab = "", xaxt = "n", yaxt = "n")
for (i in 1:11) {
  polygon(c(i - 0.5, i - 0.5, i + 0.5, i + 0.5), 
          c(-0.6, hist.self.placement[i] - 0.6, hist.self.placement[i] - 0.6, -0.6), 
          border = NA, col = "gray80")
}
abline(h = 0, lty = 3, col = "gray80")
for (i in 1:11) {
  segments(i, CATE.est.ideology[i] + CATE.SE.ideology[i] * qt(0.025, model.2.ideology$df[1]), 
           i, CATE.est.ideology[i] + CATE.SE.ideology[i] * qt(0.975, model.2.ideology$df[1]))
  points(i, CATE.est.ideology[i], pch = 19)
}
axis(1, 1:11, FALSE, lwd = 0.5)
axis(2, lwd = 0.5)
mtext(c("Liberal (left)", "Conservative (right)"), 1, 1, 
      at = c(1, 11), adj = c(0, 1))
mtext("Ideological self-placement", 1, 2.5)
mtext("Treatment effect of information on the security environment", 2, 2.5)
dev.off()

#### analysis using the binary outcome (Online Appendix E) ####
## analysis of support (Table A.2, columns 1-3)
model.1.support <- lm_robust(support.bin ~ security + pressure, data, 
                             fixed_effects = block)
round(summary(model.1.support)$coefficients[, c(1, 2, 4)], 3)

## analysis of opposition (Table A.2, columns 4-6)
model.1.oppose <- lm_robust(oppose.bin ~ security + pressure, data, 
                            fixed_effects = block)
round(summary(model.1.oppose)$coefficients[, c(1, 2, 4)], 3)

#### analysis of interaction between the treatments (Online Appendix F ####
## estimation of Model (A.1) (Table A.3)
model.A1 <- lm_robust(approve ~ security * pressure, data, 
                      fixed_effects = block)
round(summary(model.A1)$coefficients[, c(1, 2, 4)], 3)

## Figure A.2
# estimation of Model (A.1) with reversed dummy variables
model.A1.rev.1 <- lm_robust(approve ~ security * I(pressure == 0), data, 
                            fixed_effects = block)
model.A1.rev.2 <- lm_robust(approve ~ I(security == 0) * pressure, data, 
                            fixed_effects = block)

cairo_pdf("Figure_A2.pdf", width = 4, height = 3, pointsize = 7)
layout(matrix(1:2, 1, 2))
par(mar = c(3.5, 3.5, 2.5, 1), lwd = 0.5)
plot(NULL, NULL, bty = "n", xlim = c(0.5, 2.5), ylim = c(-0.4, 0.3), 
     main = "Information on\nthe security environment", 
     xlab = "", ylab = "", xaxt = "n", yaxt = "n")
abline(h = 0, lty = 1, col = "gray80")
abline(h = seq(-0.4, 0.3, 0.1)[-5], lty = 3, col = "gray80")
segments(1, model.A1$conf.low["security"], 
         1, model.A1$conf.high["security"])
points(1, model.A1$coefficients["security"], pch = 19)
segments(2, model.A1.rev.1$conf.low["security"], 
         2, model.A1.rev.1$conf.high["security"])
points(2, model.A1.rev.1$coefficients["security"], pch = 19)
axis(2, lwd = 0.5)
mtext(c("No", "Yes"), 1, 0, at = 1:2)
mtext("Information on\ninternational pressure", 1, 2.5)
mtext("Treatment effect", 2, 2.5)
plot(NULL, NULL, bty = "n", xlim = c(0.5, 2.5), ylim = c(-0.4, 0.3), 
     main = "Information on\ninternational pressure", 
     xlab = "", ylab = "", xaxt = "n", yaxt = "n")
abline(h = 0, lty = 1, col = "gray80")
abline(h = seq(-0.4, 0.3, 0.1)[-5], lty = 3, col = "gray80")
segments(1, model.A1$conf.low["pressure"], 
         1, model.A1$conf.high["pressure"])
points(1, model.A1$coefficients["pressure"], pch = 19)
segments(2, model.A1.rev.2$conf.low["pressure"], 
         2, model.A1.rev.2$conf.high["pressure"])
points(2, model.A1.rev.2$coefficients["pressure"], pch = 19)
axis(2, lwd = 0.5)
mtext(c("No", "Yes"), 1, 0, at = 1:2)
mtext("Information on\nthe security environment", 1, 2.5)
mtext("Treatment effect", 2, 2.5)
dev.off()

#### analysis of interaction using kernel smoothing (Online Appendix G) ####
## estimation using kernel smoothing
kernel.smoothing <- interflex("kernel", data, 
                              "approve", "security", 
                              "self.placement", 
                              Z = c("pressure", 
                                    paste0("block.", 1:ncol(block.dummies))), 
                              X.eval = 2:10)

## Figure A.3
cairo_pdf("Figure_A3.pdf", width = 5, height = 3, pointsize = 7)
par(mar = c(3.5, 3.5, 0.5, 0.5), lwd = 0.5)
plot(NULL, NULL, bty = "n", xlim = c(0.5, 11.5), ylim = c(-0.6, 0.4), 
     xlab = "", ylab = "", xaxt = "n", yaxt = "n")
for (i in 1:11) {
  polygon(c(i - 0.5, i - 0.5, i + 0.5, i + 0.5), 
          c(-0.6, hist.self.placement[i] - 0.6, hist.self.placement[i] - 0.6, -0.6), 
          border = NA, col = "gray80")
}
abline(h = 0, lty = 3, col = "gray80")
for (i in 1:11) {
  segments(i, kernel.smoothing$est.kernel[[1]][kernel.smoothing$est.kernel[[1]][, 1] == i, 4], 
           i, kernel.smoothing$est.kernel[[1]][kernel.smoothing$est.kernel[[1]][, 1] == i, 5])
  points(i, kernel.smoothing$est.kernel[[1]][kernel.smoothing$est.kernel[[1]][, 1] == i, 2], pch = 19)
}
axis(1, 1:11, FALSE, lwd = 0.5)
axis(2, lwd = 0.5)
mtext(c("Liberal (left)", "Conservative (right)"), 1, 1, 
      at = c(1, 11), adj = c(0, 1))
mtext("Ideological self-placement", 1, 2.5)
mtext("Treatment effect of information on the security environment", 2, 2.5)
dev.off()